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Abstract 



Millimeter-sized, spherical silicate grains abundant in chondritic me- 
teorites, which are called as chondrules, are considered to be a strong 
evidence of the melting event of the dust particles in the protoplanetary 
disk. One of the most plausible scenarios is that the chondrule precursor 
dust particles arc heated and melt in the high-velocity gas flow (shock- 
wave heating model). We developed the non-linear, time-dependent, and 
three-dimensional hydrodynamic simulation code for analyzing the dy- 
namics of molten droplets exposed to the gas flow. We conflrmed that 
our simulation results showed a good agreement in a linear regime with 
the linear solution analytically derived by Sekiya et al. (2003). We found 
that the non-linear terms in the hydrodynamical equations neglected by 
Sekiya et al. (2003) can cause the cavitation by producing negative pros- 
sure in the droplets. We discussed that the fragmentation through the 
cavitation is a new mechanism to determine the upper limit of chondrule 
sizes. We also succeeded to reproduce the fragmentation of droplets when 
the gas ram pressure is stronger than the effect of the surface tension. Fi- 
nally, we compared the deformation of droplets in the shock- wave heating 
with the measured data of chondrules and suggested the importance of 
other effects to deform droplets, for example, the rotation of droplets. We 
believe that our new code is a very powerful tool to investigate the hydro- 
dynamics of molten droplets in the framework of the shock-wave heating 
model and has many potentials to be applied to various problems. 

Keywords: meteorites, Solax System origin. Solar Nebula 
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1 Introduction 



Chondrules are millimeter-sized, once-molten, spherical-shaped grains mainly 

composed of silicate material. They are abundant in chondritic mete- 
orites, which are the majority of meteorites falling onto the Earth. They 
are considered to have formed from chondrule precursor dust particles 
about 4.56 x 10^ yr ago in the solax nebula (Amelin et al. 2002); they 
were heated and molted through flash heating events in the solar nebula 
and cooled again to solidify in a short period of time (e.g., Jones et al. 
2000 and references therein). So they must have great information on the 
early history of our solar system. Since it is naturally expected that pro- 
toplanctary disks around young stars in star forming regions have similar 
dust particles and processes, the study of chondrule formation may pro- 
vide us much information on the planetary system formation. Chondrules 
have many features: physical properties (sizes, shapes, densities, degree of 
magnetization, etc.), isotopic compositions (oxygen, nitrogen, rare gases, 
etc.), mineralogical and petrologic features (structures, crystals, degrees 
of alteration, relicts, etc.), and so forth, each of which should be a clue 
that helps us to reveal their own formation process and that of the plan- 
etary system. To reveal their formation history, many works have been 
carried out observationally, experimentally, and theoretically. 

Shock-wave heating model is considered to be one of the most plausible 
models for chondrule formation (Boss 1996, Jones et al. 2000). This model 
has been investigated theoretically by many authors (Hood & Horanyi 
1991, 1993, Ruzmaikina & Ip 1994, Tanaka et al. 1998, Hood 1998, lida et 
al. 2001, Desch & Connolly 2002, Miura et al. 2002, Ciesla & Hood 2002, 
Miura & Nakamoto 2005, 2006). One of the special features of this model 
is that chondrule precursor dust particles are exposed to the high-velocity 
rarefied gas fiow. In the gas flow, the gas frictional heating takes place to 
heat and melt the precursor dust particles (the condition to melt silicate 
dust particles was derived by lida et al. 2001). Therefore, it is naturally 
thought that the ram pressure of the gas flow affects the molten dust 
particles. Susa & Nakamoto (2002) estimated the maximum radius of the 
molten droplet above which the droplet should be destroyed by the ram 
pressure. They analyzed the non-linear phenomena like a fragmentation 
by an order of magnitude estimation. On the contrary, Sekiya et al. 
(2003) derived a linear solution of the deformation and the internal flow 
of the molten droplet exposed to the rarefied gas flow. They analytically 
solved the hydrodynamical equations assuming that the non-linear terms 
of the hydrodynamical equations as well as the surface deformation are 
sufficiently small so that linearized equations are api)ropriatc. 

The hydrodynamics of the molten droplet exposed to the rarefied gas 
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flow is an important aspect of the shock-wave heating model. The inter- 
nal flow causes the homogenization of the chemical/isotopic compositions 
and temperature distribution in the molten precursor dust particles. The 
deformation of the molten droplet might result into the deformed chon- 
drules. Moreover, the fragmentation of the molten droplet gives an upper 
limit of chondrules as pointed by Susa & Nakamoto (2002). However, 
the analysis of Susa & Nakamoto (2002) did not take into account the 
hydrodynamics of the droplet and the formulation of Sekiya et al. (2003) 
cannot be applied to such non-linear phenomena like the fragmentation. 
Additionally, the effects of non-linear terms in the hydrodynamical equa- 
tions, which was neglected in Sekiya et al. (2003), might be important in 
the hydrodynamics. In order to solve the hydrodynamical equations in- 
cluding the non-linear terms, the numerical simulations can be a powerful 
tool. The purpose of this study is to develop a numerical model of the 
molten droplet exposed to the high-velocity rarefled gas flow. 

We describe the physical model in §2 and the basic equations in §3. 
The numerical model is summarized in §4. We show the calculation results 
in §5 and some discussions in §6. Finally, we make conclusions in §7. Some 
numerical techniques we use in our model are summarized in appendixes. 

2 Physical Model 

We are considering the molten silicate dust particles exposed to the high- 
velocity rarefied gas flow. These dust particles are, of course, initially 
solid and then melted due to the gas frictional heating. In this paper, we 
do not consider how the solid particles melt in the gas flow. We assume 
the completely molten droplets (no solid lumps inside) and the physical 
properties of the droplet (the coefficient of viscosity and the fluid surface 
tension coefficient) are constant. The droplet behaves as an incompressible 
fluid because the sound speed (~ a few kms~^, Murase & McBirney 1973) 
is much larger than the fluid velocity expected in the droplet (~ a few 
tens cms~^, Sekiya et al. 2003). 

The ram pressure of the gas flow is acting on the droplet surface ex- 
posed to the high-velocity gas flow. It should be noted that the gas flow 
we consider here does not follow the hydrodynamical equations in the spar 
tial scale of chondrules because the nebular gas is too rarefied. The mean 
free path of the nebula gas can be estimated by I = l/{ns), where s is 
the coUisional cross section of gas molecules and n is the number density 
of the nebular gas. In the minimum mass solar nebula model (Hayashi et 
al. 1985), the number density of the nebula gas at 1 AU from the central 
star is about n ~ 10^''cm~^. The coUisional cross section of the hydro- 
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gen molecule is roughly estimated as s ~ 10~^^ cm'^ (e.g., HoUenbach & 
McKee 1979). So, we obtain I ~ 100cm. On the contrary, the typical size 
of chondrules is about a few 100 pim (e.g., Rubin & Keil 1984). Namely, 
the object which disturbs the gas flow is much smaller than the mean free 
path of the gas. This condition corresponds to a free molecular flow in 
which gas molecules scattered and reemitted from the droplet surface do 
not disturb the free stream velocity distribution. Therefore, in our model, 
the ram pressure acting on the droplet surface per unit area is explicitly 
given by the momentum flux of the molecular gas flow — PgU^oi) where 
Pg is the nebula gas density and Vj-^i is the relative velocity between the gas 
and the mass center of the droplet. This force causes the deceleration of 
the center of mass. In the coordinate system co-moving with the center of 
mass (we adopt this coordinate in our model), the apparent gravitational 
acceleration should be considered. 

Strictly speaking, the relative velocity of a gas molecule and a surface 
fluid element is not Vrai, since the latter has a non-zero velocity in the 
rest frame of the mass center of the droplet. However, the liquid velocity 
is much smaller than the gas molecular velocity, and is negligible (Sekiya 
et al. 2003). In hypersonic free molecular flow, the thermal and reflected 
velocities of gas molecules are also sufficiently small and are ignored. 

3 Basic Equations 

3.1 Equation of Continuity 

Generally, Eulerian methods like the CIP scheme (e.g., Yabe et al. 2001) 
use a color function (j) in the multi-phase analysis. For example, when a 
fluid of phase k occupies a certain region, the color function takes (jjk ~ 1 
inside the region and (fjk = outside. In this paper, we define tf) — 1 for 
inside the molten dust particle and (j) — for outside (ambient region). 
Using the color function, the density of the fluid element p is given by 



where the subscripts "d" and "a" indicate the inherent values for the 
molten silicate dust particle (droplet) and the ambient region, respectively. 
The other physical values of the fluid element (the coefBcient of viscosity 
fi and the sound speed Cs) are given in the same manner, namely, u = 



/id X <j!) + Pa X (1 — (^) and Cs = Cg^d x 4> + Cs,a x (1 — (/>), respectiveljl^. 

^In the numerical simulation, an undesirable situation that the hydrostatic pressure p 
becomes much larger than the value of the ambient region in spite of (/> <^ 1 could occur. 
When it occurred, we set the sound speed of the fluid element Cg — > Cg d iu order to avoid the 
difficulty in the numerical simulation. 



p = Pd X (/> + Pa X (1 - (^), 



(1) 
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ir = 0. (4) 



We need an equation that describes the time evolution of the color 
function 0. Substituting Eq. ([1} to the equation of continuity, 

^ + V-{pu)=0, (2) 

where u is the velocity, we obtain 

^ + V-{<Pu)^--P^V-u. (3) 

Integrating Eq. (|3]) over the closed region with the boundary condition of 
u = 0, we obtain 

d_ 

dt 

Therefore, it is found that the total amount of the color function 4> in the 
closed region is conserved. Then, we rewrite Eq. ([Sj as 

g + (w.V)<^ = -fl + ^l-iVv.«. (5) 

at \ Pi - 4) J 

The typical density of the silicate melt is pd ~ 3gcm^'^, while the density 
of the ambient region is extremely low (e.g., typical nebular gas density is 
about 10~^ gcm~^). Therefore, we obtain Pa/(pd — Pa) <C 1. Additionally, 
we can consider that the value of the color function (j) is order of unity 
because we especially focus on the hydrodynamics of the molten silicate 
dust particle in this study. Therefore, the second term on the right hand 
side of Eq. ((5} can be ignored. Finally, we obtain the different form of 
the equation of continuity as 

^ + {u-V)(l>^ -cpV ■ u. (6) 

It should be noted that the total amount of (p is conserved even in the 
approximated form. Eq. ([Sjl gives the time evolution of the color function 
(j). After calculating (j>, we can obtain the density p by Eq. ((T}. 



3.2 Equation of Motion 

The local velocity of the fluid element is changed by the pressure gradient, 
the viscous force, the surface tension, and the ram pressure of the high- 
velocity molecular gas flow. The ram pressure exerted on the surface of 
the droplet is given by (Sekiya et al. 2003) 

Fg = — pfm(ni ■ ng)ng5{r — n) for m ■ < 0, (7) 

where ni is the unit normal vector of the surface of the droplet, rig is the 
unit vector pointing the direction in which the gas flows, and is the 
position of the liquid-gas interface. The delta function S(r — r,) means 
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that the ram pressure works only at the interface. The ram pressure does 
not work for ni • ng > because it indicates the opposite surface which 
does not face the molecular flow. The ram pressure causes the deceleration 
of the center of mass of the droplet. In our coordinate system co-moving 
with the center of mass, the apparent gravitational acceleration g should 
appear in the equation of motion. 

The surface tension is given as (e.g., Brackbill et al. 1992) 

= -7Kni5(r-ri), (8) 

where 7 is the fluid surface tension coefficient and k is the local surface 
curvature. 

Finally, the equation of motion is written in a form 

— + {u-V)u = — ^ ^+g, (9) 

ot p 

where p is the pressure and /i is the coefficient of viscosity. In deriving the 
viscous term ^/S.u, we assume that the viscous coefficient /i do not change 
noticeably and the contribution of V ■ u in the viscosity is negligibly small 
in the droplet (Landau & Lifshitz 1987). 



3.3 Equation of State 

We can obtain the equation that describes the time evolution of the pres- 
sure p from the equation of state, which is given by 

^=c., (10) 

where Cs is the sound speed. We can rewrite it into 

dp 2 dp 

—r = Cg — , (11) 

where d/dt — d/dt+{u ■ V). Substituting the mass conservation equation 
(Eq. [2} , we obtain 

^ + {u-V)p = -pc!v ■ u. (12) 



4 Numerical Model 
4.1 Droplet and Gas Flow 

We consider the molten silicate dust particle exposed to the high-velocity 
rarefied gas flow. As we mentioned in ^ the gas flow is the free molecular 
flow because the mean free path of the gas molecules is larger than the 
typical particle size. The gas flow comes from the upstream and terminates 
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when it collide with the dust surface. Figure[T]shows the schematic picture 
of our numerical model and the coordinate system. The black thick arrows 
show the stream line of the gas flow. The gas flow gives the momentum 
to the dust particle at the area where it terminates. Therefore, the ram 
pressure exerted on the dust surface can be calculated easily, even if 
we do not solve the hydrodynamic equations for the ambient region. The 
numerical treatment of the ram pressure is described in qA.2.21 

Of course, we know that the gas flow has not only the free molecular 
motion (see Fig. [TJ, but also the random velocity caused by the thermal 
motion. It causes the hydrostatic pressure at the dust surface. Ifowever, 
the hydrostatic pressure of the ambient region is so small that it can be 
disregarded to the gas ram pressure. The typical value for the hydro- 
static pressure in the solar nebula is about psN = nksT ~ 4dynecm~^ 
for n = 10^* cm~^ and T — 300 K, where n is the number density of the 
nebula gas and T is the gas temperature. On the contrary, the ram pres- 
sure of the gas flow is about pfm ~ 4000 dyne cm~^ for the shock- wave 
heating model for chondrule formation (e.g., Uesugi et al. 2003). Addi- 
tionally, the hydrostatic pressure of the nebula gas is also much smaller 
than that inside the molten droplet. We have the well known equation for 
the hydrostatic pressure inside the droplet as po = 27/ro, where ro is the 
droplet radius. Substituting 7 — 400 dyne cm~^ for the molten silicates 
(Murase & McBirney 1973) and ro = 500 /im for the typical chondrule 
radius, we obtain po = 1.6 x lO^dynecm"^. Therefore, the hydrostatic 
pressure of the nebula gas is much smaller than that inside the molten 
droplet and the ram pressure of the gas flow. Therefore, we consider that 
the effect of psN is negligibly small for the hydrodynamics of the molten 
droplet and neglect that. 
[Figure [l] 

4.2 Ambient Region 

As we mentioned in the previous subsection, what we should do for inves- 
tigating the dynamics of the molten droplet exposed to the free molecular 
flow is just to solve the hydrodynamical equations for the molten droplet. 
Namely, there is no need to solve the hydrodynamical equations for the 
ambient region because the gas ram pressure exerted at the droplet surface 
can be calculated without doing that. 

In that case, how can we treat the outside of the droplet? In our 
numerical model, we neglect the influence of the ambient region on the 
dynamics of the molten droplet (see i|4.1|) . The best way to do that is 
to place nothing in the ambient region. It means that the fluid density 
outside the droplet is zero (see Eq. [l}. However, the equation of motion 
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cannot be solved if p = (see Eq. [Qj. Therefore, we must put the 
density of the ambient region pa > in order to solve the hydrodynamical 
equations without any other special techniques. At the same time, pa 
should be small enough not to affect the dynamics of the molten droplet. 

Additionally, we should pay attention to the sound speed of the am- 
bient region, Cs,s,. If Cs,a is large enough for the ambient region to behave 
as an incompressible fluid, the influences of the boundaries of the com- 
putational domain would travel to the molten droplet in an instant. We 
cannot say that this situation correctly simulates the molten droplet in the 
shock-wave heating. On the contrary, if Cs,a is small enough to compress 
(or expand) the ambient region without any restitution, the dynamics of 
the droplet are not affected by the boundaries. Therefore, in our numerical 
model, Cs,a should be small. 

We put Pa — 10~^gcm~^ and Ca,a — 10~^cms~^ as standard values 
for the numerical simulations (see iiA.4p . 

4.3 Coordinate system 

We set the molten droplet at the center of a cubic computational domain 
(see Fig. [!}. We adopt the Cartesian coordinate system {x,y,z), which 
is co-moving with the mass center of the droplet. Since the ram pressure 
of the gas flow explicitly appears in the equation of motion (Eq. [9l , the 
ambient region (0 = 0) does not play any important role to the droplet. 

In our model, since we examine the hydrodynamics of the droplet, it 
is no use to calculate outside of the droplet. It indicates that we con- 
sume some amount of the computational time to calculate the ambient 
region which is not important in our purpose. However, above model 
allows us to adopt the unified numerical procedure for compressible and 
incompressible fluid, which has been developed based on the concept of 
the CIP scheme (Yabe & Wang 1991). Therefore, we choose to sacrifice 
some computational time in order to make the coding process simple. 

4.4 Scheme 

In the numerical scheme, we especially should pay attention to follow- 
ing points. First, the numerical scheme for the equation of continuity 
should guarantee the mass conservation. If not, the droplet radius might 
be changed by some numerical effects. Second, the discontinuity in the 
profile of (j) between the droplet and the ambient region should be kept 
as the sharp profile. If some numerical effects make the interface diffuse, 
it would be difficult to introduce the surface tension. Final point is the 
incompressiblity of the droplet. 



10 



The equation of continuity in the form of Eq. (|6]) can be rewritten as 

1^ + V ■ icbu) = 0, (13) 

where the total amount of in a closed region should be conserved. The 
R-CIP-CSL2 scheme is one of the recent versions of the CIP scheme, 
which guarantees the exact conservation even in the framework of a semi- 
Lagrangian scheme (Nakamura et al. 2001). In order to solve Eq. (I13p . 
we adopt the R-CIP-CSL2 scheme and briefly describe this scheme in 
apDendix lA.1.21 Additionally, we make a proper correction for the velocity 
u before solving Eq. (|13|) in order to keep the incompressibility of the 
droplet (see appendix IA.3|I . 

The other two hydrodynamical equations can be separated into the 
advection phase and the non-advection phase. The advection phases of 
the equation of motion (Eq. [9| and the equation of state (Eq. I12|) are 
written as 

^ + («-V)« = 0, g + («.V)p = 0. (14) 

We solve above equations using the R-CIP scheme, which is the oscillation 
preventing scheme for advection equation (Xiao et al. 1996b, see appendix 
lA.l.ip . The non-advection phases can be written as 

du Vp Q dp nK\ 

777 = 1 . 77r=-pc,V-u, (15) 

at p p at 

where Q is the summation of forces except for the pressure gradient. 
The problem intrinsic in incompressible fluid is in the high sound speed 
in the pressure equation. Yabe & Wang (1991) introduced an excellent 
approach to avoid the problem (see appendix IA.2.H) . It is called as the 
C-CUP scheme (Yabe et al. 2001). 

Additionally, we introduce the anti-diffusion technique when calculat- 
ing the equation of continuity by the R-CIP-CSL2 scheme. The original 
R-CIP-CSL2 scheme provides an excellent solution for the conservative 
equation (e.g., Nakamura et al. 2001). However, if the initial profile 
of (j) has a sharp discontinuity (e.g., rectangle wave), the profile is los- 
ing its sharpness as the time step progresses. It is a result of the nu- 
merical diffusion. In order to keep the discontinuity of the profile, the 
anti-diffusion modification is an useful technique (e.g., Xiao & Ikebata 
2003). In this paper, we explicitly add the diffusion term with a negative 
diffusion coefficient to the CIP-CSL2 scheme (see appendix IA.1.3p . We 
perform the one-dimensional conservative equation with or without the 
anti-diffusion technique in order to show the validity of this method (see 
appendix IA.1.4|) . 

For the multi-dimensions, we use a directional splitting technique to 
perform sequential one-dimensional procedure in each direction. We adopt 
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the staggered grid for digitizing the physical variables in the hydrody- 
namical equations (see Fig. [Sjl. The scalar variables are defined at the 
cell-center and the velocity u — {u,v,w) is defined on the cell-edge. The 
surface tension and the ram pressure are defined at the cell-center. The 
numerical model of these two forces are shown in appendixes IA.2.21 and 
IA.2.31 When calculating these two forces, the gradient field of (p is re- 
quired and needs to be artificially smoothed (e.g., Yabe et al. 2001, see 
appendix IA.2.41 in this paper) . Note that the smoothing is done only 
at the calculation of these two forces and the original profile with sharp 
discontinuity is unchanged. 
[Figure [2] 

5 Results 

5.1 Input Parameters and Initial Settings 

We investigate various cases about the droplet radius, ro = 100, 200, 
500, 1000, 2000, 5000 /im, 1cm, and 2 cm. The momentum fiux of the 
gas fiow is set as pfm = 4000 dyne cm^^, which may be realized in the 
shock- wave heating model for relatively high-density shock waves (Uesugi 
et al. 2003). The coefficient of viscosity of silicate melts strongly 
depends on the temperature and the chemical composition. We adopt 
fid = 1.3gcm~^ s~^ from the model of Uesugi et al. (2003), in which they 
calculated the viscosity by the formulation of Bottinga & Weill (1972) 
assuming the temperature ~ 1800 K and the chemical composition of BO 
type chondrule. The surface tension 7 — 400 dyne cm^^ and the sound 
speed of the molten droplet Ca,d = 2 x lO^cms"^ are the typical value 
of the silicate melts (Murase & McBirney 1973). Other input parameters 
are listed in Table [1] When we consider the fragmentation process of 
the droplets in the gas fiow, it is useful to introduce the non-dimensional 
parameter We = piuiro/^, which is called as the Weber number and indi- 
cates the ratio of the ram pressure of the gas flow to the surface tension 
of the droplet. Some experiments suggested that the droplets fragment 
for We ~ 6 - 22 or higher (Bronshten 1983). 
[Table [1] 

We assume that the gas flow suddenly affects the initially spherical 
droplet (see the panel (a) in Fig. [3]|. The horizontal and vertical axes 
are the x- and y-axes normalized by the initial droplet radius ro. The gas 
flow comes from the left side of the panel. The color contour indicates the 
hydrostatic pressure p in the unit of dynecm"^ and the arrows are the 
velocity field. We have the well known equation for the initial pressure 
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inside the droplet as p — 2y/ro and p — for outside. The red curves 
indicate the contour of the color function </!> after smoothing. Thick solid 
curve indicates (j) ~ 0.5, which means the interface between the droplet 
and the ambient region, and the dashed and the dotted-dashed curves are 
(l> — 0.1 and 0.9, which are drawn in order to show the effective width of 
the transition region between the molten droplet and the ambient region. 
The initial internal velocity of the droplet is set to be all zero. The 
computational grid points are set to be 60 x 60 x 60 for the standard 
calculations. The boundary conditions are u — 0, p = 0, p = pa, </> = 0, 
and a = 0. 

5.2 Time Evolution 

Figure|3]shows the time evolution of the molten droplet exposed to the gas 
flow. The initial droplet radius ro is set as 500 jum. Since this condition 
corresponds to the Weber number We — 0.5, the droplet does not undergo 
the fragmentation (see H5.ip . A panel (a) shows the initial condition for 
the calculation. The gas flow is coming from the left side to the right side. 
The panel (b) shows the results after 0.3 x 10""^ sec. The left side of the 
droplet is directly facing the gas flow, so the hydrostatic pressure at the left 
side becomes higher than that of the opposite side. The fluid elements at 
the surface layer are blown to the downstream, but the inner velocity turns 
to upstream of the gas flow because the apparent gravitational acceleration 
takes place in our coordinate system. The droplet continues to deform 
more and more, and after 1.0 x 10~^ sec, the magnitude of the deformation 
becomes maximum. After that, the droplet begins to recover its shape to 
the sphere due to the surface tension. The panel (e) shows the shape in 
the middle of returning to a sphere, and the recovery motion is all but 
almost over at the panel (f). The droplet would repeat the deformation 
by the ram pressure of the gas flow and the recovery motion by the surface 
tension until the viscosity dissipates the internal motion of the droplet. 
[Figure [3] 

5.3 Deformation at Steady State 

As shown in Fig. [31 the droplet shows the vibrational motion: deformation 
by the ram pressure of the gas flow and recovery motion by the surface 
tension. It is expected that the vibrational motion ceases by the viscous 
dissipation and finally settles in a steady state. We plot the time evolution 
of the axial ratio C/B, where C and B are the minor axis and middle- 
length axis of the droplet, respectively. In the case of Fig. [S] it can be 
approximately considered that C and B are the droplet radii along the 
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X- and y-Eixes, respectively (see i^6.3|l . Therefore, C/B is an indicator of 
the droplet deformation, for example, C/B = 1 approximately indicates a 
sphere and smaller C/B indicates a larger deformation. 

Figure 13] shows the time evolutions of the axial ratios C/B for various 
radii of the droplets. The horizontal axis is the time after the ram pressure 
begins to act on the droplet. The solid curves are the computational 
results and the dashed lines indicate the time-averaged values defined as 

(C/B) = ^^^j^, (16) 

where the interval of the integration over t is from the left edge to the 
right edge of each dashed line. Although some time evolutions in Fig. [4] 
indicate that the vibrational motions do not terminate at the end of each 
calculation, the time-averaged values seem to well represent the axial ratio 
C/B at the final steady states. Then, we compare these time-averaged 
ajcial ratios with the linear solution derived by Sekiya et al. (2003), which 
describes the deformation of the droplet at the steady state. In order to 
compare the linear solution with our calculations in the same manner, we 
need to obtain the moment of inertia for the linear solution. In appendix 
IB] we briefly summarize how to calculate it. 
[Figure |4l 

Figure[5]shows the comparisons of the time-averaged axial ratios (filled 
circles) with the linear solution (solid curve). The horizontal axis is the 
initial droplet radius ro and the vertical axis is the axial ratio C/B. Under 
the calculation conditions we adopted here, the linear solution cannot be 
applied for ro ^ 358 /im because the Reynolds number Re exceeds unity 
(Sekiya et al. 2003). The dashed curve is a simple extrapolation of the 
linear solution. It is found that the time-averaged values (C/B) seem to 
show a good agreement not only at the radius range in which the linear 
solution can be applied but also out of the range (ro ^ 1000 /im). If the ra- 
dius exceeds about 1000 /im, the non-linear terms in the hydrodynamical 
equations appear in the simulation results. The result of ro = 2000 /im, 
in which the vibrational motion almost terminates and the droplet shape 
settles in the steady state, shows the slightly different value of C/B com- 
paring with the extrapolation of the linear solution. The panel (f) of Fig. 
[4] shows that the droplet of ro = 5000 /^m once deforms significantly at 
2 X 10~^ sec f 6 X 10~^ sec, then the surface tension has the shape re- 
cover. However, the case of the panel (f) does not undergo the vibrational 
motion around the time-averaged value (C/B), unlike other cases. After 
the first recovery motion, the droplet shape does not deform largely again. 
At the end of this calculation, the droplet settles around C/B ~ 0.9 with 
a slightly large scatter. 
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[Figure [5] 



5.4 Possibility of Cavitation 

In Fig. O we found tliat tiie axial ratio C/B of the droplet of ro = 5000 /xm 
departs from the extrapolation of the linear solution derived by Sekiya et 
al. (2003), in which the non- linear terms of the hydrodynamical equations 
as well as the surface deformation were neglected. Here, we show the 
effects of the non-linear terms in hydrodynamics of the droplet. 

Figure |6] shows the time evolution of the droplet for ro — 5000 fim. It 
is found that the droplet deforms significantly around t = 2.8 x 10~^ sec 
(panel b) . The large amount of the fluid spreads out to the perpendicular 
direction of the gas flow (y- and z-directions). Although the center of 
the flatten droplet becomes very thin, it does not fragment directly. The 
fluid blown out at once returns to the back of the droplet and gathers at 
the center (panels c and d). The center of the droplet at which the fluid 
elements gather becomes high pressure and the pressure gradient force 
pushes the fluid elements toward right in the panel. As a result, the back of 
the droplet pops out (panel e). The bump gradually disappears with time, 
and flnaUy, the droplet settles on almost the steady state with the internal 
flow largely circulating all around the droplet (panel f). In the phase, the 
hydrostatic pressure is high at the part which is directly facing the gas flow 
and extremely low around the center of eddies of the circulating internal 
fluid motion. The final pressure distribution is qualitatively different from 
the linear solution derived by Sekiya et al. (2003), in which the pressure 
at the front of the droplet (pointed by "A" in the panel) is lower than 
that at the side (pointed by "B"). 

Here, we would like to point out the behavior of the low-pressure region 
at the center of the eddies (pointed by "B" in the panel f). Generally, 
it is considered that the boiling (or vaporization) would take place in 
any liquids where the vapor pressure of the liquid exceeds its hydrostatic 
pressure. The vapor pressure of forsterite, which is one of the main com- 
ponents of chondrules, is about 1 dynecm"^ at 1850 K (Miura et al. 2002). 
On the contrary, since the hydrostatic pressure at the center of eddies is 
almost zero, the boiling might occur at the center of eddies. We call this 
phenomenon the cavitation. It indicates that the droplet might be cut off 
at the place where bubbles are generated by the cavitation. Therefore, 
it is thought that the droplet fragments into smaller pieces in the real 
situation. 

Miura et al. (2002) discussed that the ram pressure of the gas flow 
produces the high-pressure environment in the molten droplet, and the 
pressure is high enough to suppress the boiling in it when the gas frictional 
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heating is strong enough to meh the silicate dust particles. However, they 
did not take into account the hydrodynamics of the molten droplet. The 
fragmentation by cavitation might be a new mechanism to explain the 
maximum size of chondrules (see H6.ip like as the other mechanisms; the 
fragmentation due to the gas flow directly (Susa & Nakamoto 2002) or the 
stripping of the liquid surface in the gas flow (Kato et al. 2006, Kadono 
& Arakawa 2005). To simulate the behavior of the bubbles generated in 
the liquids should be the challenging task, but we think that such kind 
of simulations are very important. It might be related with the existence 
of holes in some chondrules (e.g., Kondo et al. 1997). It should be taken 
into consideration in the future work. 
[Figure [6] 

5.5 Fragmentation 

For the droplet with larger size in which the surface tension cannot keep 
the droplet shape against the gas ram pressure, the fragmentation will 
occur. Figure[7]shows the three-dimensional views of the break-up droplet. 
The droplet radius is ro = 2 cm, which corresponds to We ~ 20. In 
this case, we take a wider computational domain (—3 < x/ro < 3 and 
—4 < y/ro,z/ro < 4) than previous simulations in order to treat the 
break-up phenomenon. The computational grid points are 60 x 80 x 80, 
so the spatial resolution is worse than that of the previous cases. The 
gas flow comes from the left side of the view along the x-axis. The panel 
(a) shows the droplet shape just before the fragmentation. It is found 
that the droplet surface which is directly facing to the gas flow is stripped 
off backward. The panel (b) is just after the fragmentation. We found 
that the parent droplet breaks up to many smaller pieces. If these pieces 
collide again behind the parent droplet, the compound chondrules, which 
are thought to have collided two independent chondrules, might be formed. 

In Fig. [71 the fragmentation does not seem to be axis-symmetric. 
The reason is thought to be the directional splitting technique for multi- 
dimensions in the Cartesian coordinate. Therefore, we would not be able 
to quantitatively discuss the sizes and number of the fragments. In order 
to do that in the Cartesian coordinate system, we might need higher 
spatial resolutions. 

[Figure \7\ 
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6 Discussions 



6.1 Maximum Size of Chondrules 

Susa & Nakamoto (2002) suggested that the fragmentation of the droplets 
in high- velocity gas flow limits the sizes of chondrules (upper limit). They 
considered the balance between the surface tension and the inhomogene- 
ity of the ram pressure acting on the droplet surface, and derived the 
maximum size of molten silicate dust particles above which the droplet 
would be destroyed by the ram pressure of the gas flow using an order of 
magnitude estimation. In their estimation, they adopted the experimen- 
tal data in which the droplets suddenly exposed to the gas flow fragment 
for We <; 6 (Bronshten 1983, p. 96). This results into the fragmenta- 
tion of droplet for ro 6000 /im if we adopt our calculation conditions 
(Pfm ~ 4000 dynecm"^ and 7 = 400 dynecm"^). 

On the contrary, our simulations showed the possibility that the droplets 
fragment through the cavitation even if the radius is smaller than the max- 
imum size derived by Susa & Nakamoto (2002) (see i]5.4p . The cavitation 
would occur when the pressure gradient force cannot support the fluid 
motion to circulate around the eddy. Here, we consider the condition 
that the cavitation would take place using a simple order of magnitude 
estimation. The pressure gradient force per unit volume can be estimated 
as /pres ~ p/(ro/2), where ro/2 indicates the radius of the eddy. On the 
other hand, the centrifugal force for the fluid element circulating around 
the eddy can be given as /cent ~ Pd'^circ/(^o/2), where v^rc is the rota- 
tional velocity. Balancing the pressure gradient and the centrifugal force 
around the eddy, we have 



Substituting p — 27/ro and Vd^c — ^max = 0.112pfniro//id (Sekiya et al. 
2003), we obtain 



This equation gives the critical radius of the droplet above which the 
cavitation might occur in the center of the eddy and result into the frag- 
mentation. 

Figure [8] shows the two critical radii as a function of the viscosity; 
the radius for We — 6 (dotted line) and the radius at which the pres- 
sure gradient force balances with the centrifugal force (dashed line). We 
assume Pfm = 4000 dynecm"^, 7 — 400dynecm~^, and pd ~ 3gcm~^, re- 
spectively. Our simulation results are also shown; not fragment (circles), 
cavitation would take place (flUed circles), and fragmentation (triangle). 




(18) 
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respectively. It is found that there is a critical value for the viscous co- 
efBcient of the molten droplet /icr below which the cavitation would take 
place at a certain region in radius, even if We is not so large that the 
fragmentation is not expected to occur. The critical viscous coefficient 
/icr is calculated as 

/6Mm2Vpd 

\ iPfm 

= 12.8 ( ^ ( 2 ,] gcm-i s-(l9) 

\ 4000 dyne cm- 2 y \^400dynecm"i / ^ ^ 

Susa & Nakamoto (2002) clearly showed the maximum sizes of molten 
droplets in the post-shock region, as a function of the initial gas den- 
sity po and the shock velocity Vs (see Fig. 1 in their paper). They also 
compared flcr with the shock conditions for the chondrule formation (in 
which the chondrule precursor dust particles can melt and do not evap- 
orate completely). In their figure, they assumed that the droplet will be 
destroyed by the gas ram pressure when We > 6. However, if the droplet is 
so heated that the droplet viscosity exceeds the critical value /Xcr, even the 
smaller droplet would be destroyed by the cavitation. From above spec- 
ulations, it is expected that BO chondrules (which are believed to have 
melted totally and experienced high temperature, e.g., Hewins & Radom- 
sky 1990) are smaller than POP chondrules (which are thought to have 
melted partially). On the contrary, there are no obvious size differences 
between BO and POP chondrules (e.g., Rubin 1989). This discrepancy 
might be explained by that the precursor dust particles were not so large. 
If the maximum precursor size is smaller than about 1000 /im, the cav- 
itation would not take place for jj, Igcm^^s"^ in Fig. [S] Another 
possibility is that the viscosity of the molten droplets have not been be- 
low the critical value Hcr, above which the cavitation would not take place. 
Namely, if the dust parameters do not enter the region of "cavitation" in 
Fig. [51 the size-dependence in chondrule sizes would not appear by the 
cavitation. In order to investigate above issue, we might have to consider 
not only the hydrodynamical effect but also the thermal evolutions of the 
molten droplets in order to discuss the maximum sizes. The thermal evo- 
lution of the droplets is beyond the scope of this paper. We would like to 
investigate this issue in the future work. 
[Figure [5] 

6.2 Final Shape of Droplet 

In this paper, we numerically simulated the hydrodynamics of molten 
droplets exposed to the gas ffow. We assumed that the coefficient of 
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viscosity is constant in the droplets, and the droplets melt completely 
and the physical conditions inside the droplets are homogeneous. On 
the contrary, the shapes of chondrules are considered to be determined 
around the re-solidification phase. Before the molten droplets re-solidify, 
the droplet temperature falls down and the coefficient of viscosity would 
change significantly. Moreover, if the droplets melt partially, the physical 
conditions inside the droplets should not be homogeneous. Therefore, we 
must discuss how we can consider above effects on the final shapes of 
chondrules. 

It is naturally considered that the viscosity of the droplet becomes 
larger as the droplet temperature falls down. In such a lower Reynolds 
number environment, the linear solution by Sekiya et al. (2003) is a 
better approximation to describe the external shape and internal fiow of 
the droplet. According to their linear solution, the external shape of the 
droplet does not depend on the viscosity. It indicates that the external 
shape is determined mainly by the balance between the surface tension 
and the gas ram pressure. Therefore, it is considered that the droplet 
shape does not change significantly at the cooling phase in which the 
viscosity of the droplet becomes higher and higher. 

Another problem that we must take care for considering the droplet 
shapes is inhomogeneity of the physical conditions in the droplet. For 
example, if the dust particle is not heated enough, it would not melt 
completely and the solid lumps floating in the droplet would exist. In 
this case, the assumption of the completely molten droplet we adopted 
loses the validity. However, if the external shape is significantly affected 
by the solid lumps, the surface of chondrules should be irregular and 
cannot be approximated by smoothed shapes (e.g., spheres or ellipsoids). 
If considering conversely, it is thought that we can remove the above effect 
by considering only the chondrules which have smooth surfaces. In fact, 
there are some fraction of chondrules with smooth surfaces (Tsuchiyama et 
al. 2003). Such chondrules would have melted almost completely, at least 
the unmelted solid lumps do not affect the external shapes of droplets. 

However, the majority of chondrules has the porphyritic texture, which 
indicates that most of chondrule precursors did not completely melt in 
the heating event (Hewins & Radomsky 1990). Moreover, iron sulfide 
inclusions are observed in natural chondrules with various forms (Uesugi 
et al. 2005). The unmelted clumps and the iron sulfide inclusions might 
disturb the internal ffow of the molten silicate particle^. This problem 

^Uesugi et al. (2005) considered that the trajectories of the iron sulfide inclusions in 
the silicate melts, however, they did not take into account the influence of the iron sulfide 
inclusions on the flow pattern of the ambient silicate melt. 



19 



is very complex and beyond the scope of this paper. We beheve that this 
problem can be investigated in the future. 

6.3 Comparison with Chondrules 

Tsuchiyama et al. (2003) studied three-dimensional shapes of chondrules 
using X-ray microtomography. They measured twenty chondrules with 
perfect shapes and smooth surfaces, which were selected from 47 chon- 
drules separated from the Allende meteorite (CVS). The external shapes 
were approximated as three-axial ellipsoids with a-, b-, and c-axes (axial 
radii are A, B, and C (A > B > C), respectively) using the moments of 
inertia of the chondrules, where the rotation axes with the minimum and 
maximum moments correspond to the a- and c-axes, respectively. They 
found that (1) the shapes are diverse from oblate (A ~ B > C), general 
three-axial ellipsoid (A > B > C) to prolate chondrules (A > B ~ C), 
and (2) two groups can be recognized: oblate to prolate chondrules with 
0.9 < B/A < 1.0 (group-A) and prolate chondrules with 0.74 < B/A < 
0.80 (group-B). 

The axial ratio C/B of the group-A chondrules is about 0.9 ^ C/B <, 
1.0 and the radii are about from 200 /im to 1000 /xm (A. Tsuchiyama, 
private communication). On the contrary, for the molten droplet exposed 
to the gas flow, C/B also depends on the momentum flux of the rarefied 
gas flow pfm. Figure [9] shows the axial ratio C/B for various values of 
Pfm (= 400, 1000, and 4000 dyne cm~^) calculated by the linear solution 
derived by Sekiya et al. (2003). We also display the data of the group-A 
chondrules; filled circles are oblate (B/A > C/B) and open circles are 
others. The radii of chondrules are estimated by (ABC)^^"^. It is found 
that the data of group-A chondrules seem to distribute between the linear 
solutions for = 400 dyne cm~^ and — 4000 dyne cm~^. In other 
words, assuming that the axial ratio C/B of those chondrules were caused 
by the ram pressure of the gas flow when they melted, the expected ram 
pressure of the gas flow is about 400 — 4000 dyne cm~^. This value is the 
typical value of the shock-wave heating model for chondrule formation. 
Therefore, it is strongly suggested that the group-A chondrules might 
have been formed in the high-velocity rarefied gas flow of the shock-wave 
heating. 

[Figure [9] 

However, it should be noted that some of those group-A chondrules 
have the axial ratios B/A which are not unity. This indicates that the 
shapes are not pure oblate but general three-axial ellipsoids. Moreover, 
the group-B chondrules show relatively small values of B/A (~ 0.75), 
but C/B ~ 1. These shapes are called as "prolate." Such not-oblate chon- 
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drules cannot be explained considering only the effect of the ram pressure. 
In order to produce the not-oblate molten droplet, we must investigate 
other effects which give some dynamical effect to the droplet (e.g., rotation 
of the droplet). In our forthcoming paper, we are planing to investigate 
the rapidly rotating droplets exposed to the gas flow. The dust rotation 
causes the three-dimensional deformation depending on the rotation axis. 
Therefore, the analysis based on the assumption of axis-symmetry would 
not be used. Since our numerical code has been developed on the concept 
of the three-dimension, any special improvements are not required for the 
application to the rapidly rotating droplets. 

6.4 Ram Pressure at Re-solidification? 

The precursor dust particles moving through the nebula gas are deceler- 
ated by the ram pressure of the gas flow. Therefore, the relative velocity 
Vre\ between them becomes small with time. It indicates that the gas ram 
pressure Pfm = Pg^rei a-t the time when the molten droplet re-solidifies 
would be smaller than that when it melts first. 

Figure [10] shows the dynamical/thermodynamical evolutions of pre- 
cursor dust particles in the post-shock region (Miura & Nakamoto 2006). 
The dust radius is 1000 /xm. The horizontal axis is the distance from the 
shock front x in logarithmic scale. The top panel shows the dust temper- 
ature Td and the radiation temperature defined by Trad = (t^J'/o'sb)^^'^, 
where J is the mean intensity of the radiation field mainly emitted from 
the dust particles and the shocked gas in the post-shock region, and asB is 
the Stefan-Boltzmann constant. The middle panel shows the dust veloc- 
ity Vd and the gas velocity Vg^ against the shock front. The bottom panel 
shows the ram pressure acting on the dust particles. The shock velocity 
Vb = lOkms"^, the pre-shock gas number density no = 10^'* cm~^, and 
the power-law dust size distribution are assumed. The mean intensity J 
also depends on the dust-to-gas mass ratio Cd- The left column indicates 
a case of the lower dust-to-gas mass ratio (Cd = 0.01). The dust tem- 
perature increases up to about 1700 K around x ~ 10^ km and after that 
decreases. Although the relative velocity between the dust particle and 
the gas decreases monotonically with the distance, the ram pressure pfm 
increases until x ~ 10^ km because the gas density pg increases, and after 
that decreases. However, it is found that the ram pressure does not de- 
crease significantly when the dust temperature falls down to the melting 
point of silicates (1573 K, thin dashed line, e.g., Tachibana & Huss 2006). 
On the contrary, the right column is a more dusty case (Cd ~ 0.10). In 
that case, the radiation temperature is stronger than the left panels as a 
result of the blanket effect and in this case, exceeds the melting point of 
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silicates. Since the dust particles are also heated by the radiation field, 
they do not cool even if the relative velocity is almost zero. They can 
cool if the radiation field becomes weak. At that time, the ram pressure 
has not already taken place on the dust particles. To summarize, the ram 
pressure acting on the re-solidifing molten droplet can take various values 
depending on the shock conditions and dust models in the chondrule- 
forming region. Generally, the low dust-to-gas mass environment would 
result into the large ram pressure. 
[Figure [iD] 

6.5 Dust Particle Rotation 

We investigated the hydrodynamics of molten droplet exposed to the gas 
flow. However, in the shock-wave heating model, the situation where the 
rapidly rotating droplets are exposed to the high-velocity gas flow can be 
considered. Some possibilities for the origin of the dust rotation can be 
considered, for example, the net torque when the droplet fragments by 
the gas ram pressure (Susa & Nakamoto 2002), the interaction between 
the dust particles and the ambient gas flow, and the dust-dust collision. 
Therefore, we have to take into account the effect of the rotation. When 
the rotation axis is perpendicular to the gas flow, the three-dimensional 
effect is expected in the droplet shapes. Our code has been already devel- 
oped in the concept of the three-dimensional calculations. We are planing 
to examine the case that the rotating droplets are exposed to the high- 
velocity gas ffow in our forthcoming paper. 

7 Conclusions 

We have developed the non-linear, time-dependent, and three-dimensional 
hydrodynamic simulation code in order to investigate the hydrodynamics 
of molten droplets exposed to the high-velocity rarefled gas flow. Our 
numerical code is based on the concept of the CIP scheme, guarantees the 
exact mass conservation, and additionally, we introduced the anti-diffusion 
technique for suppressing the effect of the numerical diffusion at the dis- 
continuity in density. We carried out the simulations for various droplet 
radii (ro — 100 /im — 2 cm). The other physical parameters were adopted 
typical values for the shock-wave heating condition (the momentum of 
molecular gas flow pf,,, ~ 4000 dyne cm~^) and the silicate melts (surface 
tension 7 = 400dynecm~^, coefficient of viscosity — 1.3gcm~^ s~^). 
We conclude as follows: 

1. Since we considered that the gas ram pressure suddenly affects the 
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initially spherical droplet, the droplets whose radii are smaller than 
5000 /um showed the vibrational motion (deformation by the ram 
pressure and recovery motion by the surface tension) . The vibration 
gradually dissipates by the viscosity and the droplets tend to settle 
in the steady states. 

2. The degree of the deformation at the steady state (or time-averaged 
value) showed a good agreement with the linear solution derived by 
Sekiya et al. (2003) if the droplet radii are smaller than 1000 jim. 
For larger droplets, the final shapes are far from the linear solu- 
tion because the non-linear term in the hydrodynamical equation 
dominates, which terms were ignored in the derivation of the linear 
solution. 

3. If the droplet radii are larger than 5000 fjxn, the hydrostatic pressure 
inside the eddies of the droplet internal flow becomes almost zero. 
In this region, the phase transition from liquid to vapor would occur 
and some bubbles would be generated in the droplet (cavitation). 
It is considered that the cavitation causes the fragmentation of the 
droplets. 

4. When the droplet radius is larger than 1 cm, the droplet fragments 
directly by the gas ram pressure. We found that the parent droplet 
breaks up to many smaller pieces. If these pieces collide again behind 
the parent droplet, the compound chondrules, which are thought to 
have collided two independent chondrules, might bo formed. 

5. We considered that the fragmentation through the cavitation might 
be a new mechanism to determine the maximum size of chondrules. 
Using the order of magnitude estimation, we found that the cavi- 
tation can easily occur in the low-viscosity droplets. We discussed 
the possibility that this relation can explain why chondrules in some 
carbonaceous chondrites are smaller than that in other chondrites. 

6. We compared our simulation results of the molten droplets with 
chondrules measured by Tsuchiyama et al. (2003) in the external 
shape, and found that the variety of chondrule shapes cannot be re- 
produced only by the effect of the ram pressure of the gas flow. We 
pointed out the importance of other mechanisms, e.g., the rotation 
of the droplets. 
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A Numerical Method in Hydrodynamics 

In this appendix, we briefly explain the strategies of the numerical schemes 
that we adopted in our model. The hydrodynamical equations except for 
the equation of continuity were separated into two phases; the advection 
phase f i]A.l|l and the non-advection phase f i]A.2p . The equation of conti- 
nuity was solved in the conservative form f i]A.1.2|l . We also describe the 
method to keep the incompressibility of the fiuid in i]A.3l and the model 
parameters we adopted in this study in HA.4\ 



A.l Advection Phase 
A. 1.1 CIP scheme 

The CIP scheme is one of the high-accurate numerical schemes for solv- 
ing the advection equation (Yabe et al. 2001). In one-dimension, the 
advection equation is written as 

where / is a scaler variable of the fiuid (e.g., density), u is the fluid velocity 
toward the a:-direction, and t is the time. When the velocity u is constant, 
the exact solution of Eq. (|20p is given by f{x;t) — f{x — ut;0), which 
indicates a simple translational motion of the spatial profile of / with the 
velocity u. 

Consider that at the time step n, the values of / on the computational 
grid points Xi-i, Xi, and Xi+i (/jli, f", and /i+i) are given as the filled 
circles in Figure [TT] and the updated value at x = Xi (/"^^) is now re- 
quired, where the updated value indicates the value at the next time step 
n + 1. The time interval between these time steps is set to be At. From 
the solution of Eq. (I2Q|I . we can obtain /"^^ by just calculating /" at 
the upstream point x — Xi ~ uAt. If the upstream point locates between 
Xi-i and Xi, we have to interpolate /" with an appropriate mathematical 
function composed of fP-i, and so forth. There are some variations of 
the numerical solvers by the difference of the interpolate function Fi{x). 
One of them is the first-order upwind scheme, which interpolates /" by 
a linear function and satisfies following two constraints; Fi{xi-i) = f"_i 
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and Fi{xi) = /f. The other variation is the Lax-WendrofF scheme, which 
uses a quadratic polynomial satisfying three constraints; Fi{xi^i) — f"-i, 
Fi(xi) = /r,andF,(i,+i) = /r+i. 

On the contrary, the CIP scheme interpolates values using a cubic 
polynomial, which satisfies following four constraints; Fi{xi-i) — f"-i, 
Fi{xi) = /r, dFi/dx{xi-i) = and dF/dxixi) = where = 

df/dx is the spatial gradient of /. The way to update fx is as follows. 
Differentiating Eq. (I20p . we obtain 

dfx ,dfx du 

where the second term of the left hand side indicates the advection and 
the right hand side indicates the non-advection term. The interpolate 
function for the advection of fx is given by dFi/dx. The non-advection 
term can be solved analytically by considering that du/dx is constant. 
Additionally, there is an oscillation preventing scheme in the concept of 
the CIP scheme, in which the rational function is used as the interpolate 
function. The rational function is written as (Xiao et al. 1996a) 

, , _ ai(x - Xjf + bi{x - Xjf + Ci{x - Xj) + 
^^'^ l+a,P.{x-x,) ' ^^^> 

where a^, bi, Ci, Ui, and Pi are the coefficients which are determined from 
fi^-i, fr, fx.i-i, and fx^i- The expressions of these coefficients are shown 
in Xiao et al. (1996a). This scheme is called as the R-CIP scheme. 

Fig. [11] also shows the interpolate functions with various methods; 
CIP (solid), Lax-Wendroff (dashed), and first-order upwind (dotted). In 
this figure, we assume /" = at all grid points for the CIP scheme. The 
difference in the numerical solutions of Eq. (|20|l is discussed in i)A.1.4l 
[Figure [Tl] 

A.1.2 CIP-CSL2 scheme 

The CIP-CSL2 scheme is one of the numerical schemes to solve the con- 
servative equation. In one-dimension, the conservative equation is written 

Integrating Eq. (I23p over x from Xi to Xi+i, we obtain 



9<Ti+l/2 
dt 



uf 



i+1 



= 0, (24) 



where = J^'^^ fdx. Since the physical meaning of uf in the second 

term of the left hand side is the fiux of cr per unit area and per unit time, 
the time evolution of a is determined by 

""1+1/2 — ""1+1/2 — + 1 + Ji, (25) 
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where Ji = ufdt is the transported value of cr from x < Xi to x > Xi 

within Af, and it corresponds to the grayed region in Figure 1121 where 
the filled circles indicate /" and the thick solid curve is the interpolate 
function Fi{x) for a region Xi^\ < x < Xi. We assume Ui > Q m this 
figure. 

The CIP-CSL2 scheme uses the integrated function Di {x) = Fi {x)dx 

for the interpolation when Ui > 0. The function Diix) is a cubic polyno- 
mial satisfying following four constraints; Di{xi-i) — 0, Di{xi) = (Ji-1/2, 
dDi/dx{xi-i) = Fi{xi^i) = and dDi/dx{xi) = Fi{xi) = fi. More- 

over, since Eq. (|23|) can be rewritten as the same form of Eq. (|2ip . we 
can obtain the updated value /"^^ as well as f^^^ in the CIP scheme 
(^ATlTJ- Additionally, there is an oscillation preventing scheme in the 
concept of the CIP-CSL2 scheme, in which the rational function is used 
for the function Di{x). The rational function is written as (Nakamura et 
al. 2001) 

„ , X _ aijx - Xj)''^ + bjjx - Xj)^ + Ci{x - Xj) 

where a^, bi, Ci, at, and /3i are the coefficients which are determined from 
fi-i, fi, and cri_i/2- The expressions of these coefficients are shown in 
Nakamura et al. (2001). This scheme is called as the R-CIP-CSL2 scheme. 
The example of the test calculation is shown in i]A.1.4l 
[Figure [H] 

A. 1.3 Anti-DifFusion 

In order to keep the sharp discontinuity in the profile of (f}, we explicitly 
add the diffusion term with a negative diff'usion coefficient a to the CIP- 
CSL2 scheme. In our model, we have an additional diffusion equation 
about a (see Eq. I24|l as 

da d f da\ 

m = 8-x[^d-x)- (27) 

Eq. (|27p can be separated into two equations as 

where J' indicates the flux per unit area and per unit time. Using the 
finite difference method, we obtain 

<^i+l/2 = <^i+l/2 — {Ji+l — Ji), (30) 
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J'i = --&i{'^i+l/2 — (31) 

where J' = J' /{Ax /At) is the mass flux which has the same dimen- 
sion of a and a = a/{Ax'^/At) is the dimensionless diffusion coefficient. 
The superscript * and ** indicate the time step just before and after the 
anti-diffusion. However, if we use Eq. (|3ip directly, we cannot obtain 
appropriate solutions. Therefore, we make some inventions in Eq. (|3H) to 
obtain the mass flux J'. 

In this paper, we calculate J' as 

Ji = —o-i X minmod(5'i-i, 5";, ^i+i), (32) 

where Si = CTi+1/2 — o"i_i/2- The minimum modulus function (minmod) 
is often used in the concept of the flux limiter and has a non-zero value of 
sign(a) min(|a|, |&|, |c|) only when a, b, and c have same sign. The value of 
the diffusion coefficient a is also important. Basically, we take a — —0.1 
for the anti-diffusion. Here, it should be noted that a takes the limited 
value as < (T < am, where (Jm is the initial value for inside of the droplet. 
The undershoot (cr < 0) or overshoot (cr > am) are physically incorrect 
solutions. To avoid that, we replace ai — 0.1 only when (7i_i/2 or 0-^+1/2 
are out of the appropriate range. 

We insert the anti-diffusion phase after the CIP-CSL2 scheme is com- 
pleted. We also have the anti-diffusion for other directions (y and z) in 
the same manner. The test calculations for one-dimensional conservative 
equation with and without the anti-diffusion are shown in iiA.1.41 

A. 1.4 Test Calculation 

In order to demonstrate the advantage of the CIP scheme, we carried out 
one- dimensional advection calculations with various numerical schemes. 
Figure [13] shows the spatial profiles of / of the test calculations. The 
horizontal axis is the spatial coordinate x. The initial profile is given by 
the thin solid line, which indicates a rectangle wave. We set the fiuid 
velocity u = 1, the intervals of the grid points Aa; = 1, and the time 
step for the calculation At = 0.2. These conditions give the CFL number 
V = uAt/ Ax — 0.2, which indicates that the profile of / moves 0.2 times 
the grid interval per time step. Therefore, the right side of the rectangle 
wave will reach a:: = 80 after 300 time steps and the dashed line indicates 
the exact solution. The filled circles indicate the numerical results after 
300 time steps. 

The upwind scheme does not keep the rectangle shape and the profile 
becomes smooth by the numerical diffusion (panel a) . In the Lax- Wendroff 
scheme, the numerical oscillation appears behind the real wave (panel b). 
Comparing with above two schemes, the CIP scheme seems to show better 
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solution, however, some undershoots (/ < 0) or overshoots (/ > 1) are 
observed in the numerical result (panel c). In the R-CIP scheme, although 
the faint numerical diffusion has still remained, we obtain the excellent 
solution comparing with the exact solution (panel d). 

We also show the numerical results of the one-dimensional conservative 
equation (Eq. [23]). We use the same conditions with the one-dimensional 
advection equation (Eq. I20|) . Note that Eq. (I23|l corresponds to Eq. (|20[) 
when the velocity u is constant. The panel (e) shows the results of the 
R-CIP-CSL2 scheme, which is similar to that of the R-CIP scheme. In 
the panel (f ) , we found that the combination of the R-CIP-CSL2 scheme 
and the anti-diffusion technique ( i]A.1.3|) shows the excellent solution in 
which the numerical diffusion is prevented effectively. 
[Figure [13] 



A. 2 Non- Advection Phase 
A.2.1 C-CUP scheme 

Using the finite difference method to Eq. p5|l . we obtain 

X7 = — + —' — X7 — = -pc,V-u , (33) 

At p* p* At 

where the superscripts * and ** indicate the times before and after calcu- 
lating the non-advection phase, respectively. Since the sound speed can 
be very large in the incompressible fluid, the term related to the pressure 
should be solved implicitly. In order to obtain the implicit equation for 
p** , we take the divergence of the left equation and substitute u** into 
the right equation. Then we obtain an equation 

(34) 

The problem to solve Eq. p4p resolves itself into to solve a set of linear 
algebraic equations in which the coefficients become an asymmetric sparse 
matrix. After p** is solved, we can calculate u** by solving the left 
equation in Eq. (|33lFI . 




A. 2. 2 Ram Pressure of Gas Flow 

Consider that the molecular gas flows for the positive direction of the 
a:-axis. The x-component of the ram pressure Fg^x is given by 

Fg,x = PiinS(x - Xi), (35) 

•^In our model, we neglect the viscous term in Q when calculating Eq. 1341 In the original 
C-CUP scheme developed by Yabe & Wang (1991), the other terms (correspond to the surface 
tension and the ram pressure, in our model) are also ignored when calculating the pressure. 
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where Xi is the position of the droplet surface. This equation can be 
separated into two equations as 

-g^ = ~PfmO(x — Xi), (36) 

F.,. (37) 
where M is the momentum flux of the molecular gas flow. Eq. (|36p means 
that the momentum flux terminates at the droplet surface. Eq. (|37p 
means that the decrease of the momentum flux per unit length corresponds 
to the ram pressure per unit area. 

Using the finite difference method to Eq. (|36|) . we obtain 



Mi+l = Mi -P{in(4>i + 1 - 4>i) for 't>i + i > 4>i, (38) 

where cj) is the smoothed profile of (j> (see appendix lA.2.4|l , and Mi+i = Mt 
for (t>i+i < 4>i because the momentum flux does not increase when the 
molecular flow goes outward from inside of the droplet. Similarly, we 
obtain ^ 

Fg,xi = (39) 

from Eq. (|37|l . The momentum flux at upstream is Mo = Pfm- First, 
we solve Eq. (|38p and obtain the spatial distribution of the molecular 
gas flow in all computational domain. After that, we calculate the ram 
pressure by Eq. (|39p . 



A. 2. 3 Surface Tension 

The surface tension is the normal force per unit interfacial area. Brackbill 
et al. (1992) introduced a method to treat the surface tension as a volume 
force by re-placing the discontinuous interface to the transition region 
which has some width. According to them, the surface tension is expressed 
as 

= T'^Vf^/M, (40) 

where [(j)] is the jump in color function at the interface between the droplet 
and the ambient gas (in our definition, we obtain [0] = 1). The curvature 
is given by 

fc = -(V-n), (41) 

where 

n = V(^/!V0!. (42) 

The finite difference method of Eq. (|42p is shown in Brackbill et al. (1992) 
in detail. When we calculate the surface tension, we use the smoothed 
profile of 4> (see appendix lA. 2. 4p . 
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A. 2. 4 Smoothing 

We can obtain the numerical results keeping the sharp interface between 
the droplet and the ambient region. However, the smooth interface is 
suitable for calculating the smooth surface tension. We use the smoothed 
profile of (j) only at the time to calculate the surface tension and the ram 
pressure acting on the droplet surface. In this study, the smoothed color 
function (j) is calculated by 

7 1, 1 <^«.j.fc + Ci -^-c-i + C2 X^ia <7^i2 + Ca <?^i3 
'^^2'^-'^+2 1 + 6C71 + 12C. + 8C3 ' ^'^^ 

where Li, 1/2, and L3 indicate grid indexes of the nearest, second nearest, 
and third nearest from the grid point (i, j, fc), for example, Li = j, k), 
Z/2 = + j' + l, k), Z/3 = + j + fc + 1), and so forth. It is easily found 
that in the three-dimensional Cartesian coordinate system, there are six 
for Li, twelve for L2, and eight for La, respectively. The coefficients are 
set as 

Ci = 1/(6 + + 8/\/3), C2 = C1/V2, Ca = Ci/\/3. (44) 

We iterate the smoothing five times. Then, we obtain the smooth transi- 
tion region of about twice grid interval width. We use the smooth profile 
of (f) only when calculating the surface tension and the ram pressure. It 
should be noted that the original profile (j) with the sharp interface is kept 
unchanged. 

A. 3 Pressure Correction for Incompressible Fluid 

The molten chondrule precursor dust particles can be regarded as incom- 
pressible because of the large sound velocity (see ^IJ. However, in the 
numerical solutions, it can occur that the velocity in the droplet calcu- 
lated by Eq. (|33p has some divergence. In order to vanish it, we need 
to perform an extra iteration step, for example, as proposed by Chorin 
(1968). 

Although we obtained the velocity u** after calculating Eq. (I33|) . the 
divergence of u** inside the droplet is not small enough to be expected 
from its large sound speed (incompressibility) . In order to obtain the 
proper velocity, we need to correct the pressure by some method and re- 
calculate the velocity by using the corrected pressure. Here, we write the 
velocity and pressure after the correction as u** and p** . Additionally, 
we define the correction of the pressure as 5p = p** — p" . The velocity 
should be changed according to the equation of motion. On the other 
hand, the pressure should be changed according to the equation of state. 
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Therefore, we have the velocity /pressure correction method as (c.f., Eq. 
|33]) 

w** — w** V(5p p** ^ p** „ 

^^^--P^?V.. . (45) 

We iterate these equations until p** converges. The same numerical 
scheme as applied for Eq. P5)l can be used to solve Eq. ifiS)). 



A. 4 Model Parameters 

In Table [T] parameters marked by (*) have no physical meanings. These 
values are empirical settings for the numerical simulation. 

While we set the density of the ambient region pa = 10~® gcm~"^, the 
typical gas density of the protoplanetary disk is ~ lO^^gcm"'^ in the 
minimum mass solar nebula (Hayashi et al. 1985). There is a difference 
of about four orders of magnitude in the densities, however, the numerical 
results do not depend significantly on the ambient density because it is 
too small to afi'ect the dynamics of the molten droplet. The stronger 
the density contrast between the droplet and the ambient region is, the 
more difficult the numerical simulation becomes. From above reasons, we 
adopted the relatively large value for the ambient density in the numerical 
simulation. 

While we set the sound speed of the ambient region Cs,a = 10~^ cms~^, 
the typical value of the disk gas is Cg.diak — (fcBT/w-Ha)^''^ — lO^cms"^, 
where fee is the Stefan-Boltzmann constant, T is the gas temperature (we 
substitute T = 300 K), and is the molecular weight of H[2. Remember 
that in this study, the disk gas around the molten droplet behaves like as 
a free molecular fiow and does not follow the hydrodynamical equations 
(^. Therefore, we do not need to consider the sound wave and the 
pressure change by the compression/expansion in the ambient region. We 
thought that the extremely low sound speed for the ambient region is one 
of the appropriate ways to express above features in the concept of the 
multi-fluids analysis. 

Finally, we set the viscous coefficient of the ambient region = 
10~^ g cm~^ s"'^ because it should not be zero for the stability of the cal- 
culations. 



B Moment of Inertia for Linear Solution 

The external shapes obtained by our numerical simulations were approxi- 
mated as three-axial ellipsoids using the moments of inertia. We compared 
our results with the linear solution derived by Sekiya et al. (2003) in this 
paper. When we obtain the axial ratio of the linear solution, the moment 
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of inertia for the linear solution is needed. It can be calculated using the 
formalism by Sekiya et al. (2003) as follows. 

Sekiya et al. (2003) gave the distance between the dust surface to 
the origin of the spherical coordinate^ as rs{0) — ro(l + WeXsiO)), where 
Ws — pfmTo/'y is the Weber number, 9 is the angle from the radius of 
droplet to the negative x-direction in our coordinate, and Xb is the dimen- 
sionless parameter which was given as a function of 6. The expression of 
Xs is given in Sekiya et al. (2003). The moment of inertia around the 
a;-axis is calculated as 

Ixx ~ Pd dr d9 dtp {r sin 9)^ sin 6, 

Jo Jo Jo 

5 

= ^Trpdro' J2 ^Prnbn^Wr, (46) 

m = 

where P is the permutation defined by „Pr = n\/{n — r)\ and 

hm= [ xT{9)sin^9d9. (47) 
Jo 

We obtain bo — 4/3. It can be easily confirmed that I^^ has the same value 
of a sphere (/ = (8/15)7rpd'"o) when We = 0. For other coefficients, we 
calculate using an appropriate numerical integration method and obtain 
6i = 0.42222 X 10"\ 62 = 0.51562 x 10"^ 63 = 0.74446 x 10"", 64 = 
0.44106 X 10"", and 65 = -0.28485 x 10"^ 



"Sekiya et al. (2003) chose the origin of the spherical coordinates as dra/dO = at 6 = 
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Table 1: Input physical parameters for simulations of molten droplet exposed 
to the gas flow. The meanings of the parameters marked by (*) are shown in 
Appendix I A. 41 



parameter 


sign 


value 


momentum of gas flow 


Pirn 


4000 dynecm"^ 


surface tension 


7 


400 dyne cm~^ 


viscosity of droplet 




1.3gcm~^ s~^ 


density of droplet 


Pd 


3 gcm^"^ 


sound speed of droplet 




2 X lO^cms-i 


(*) density of ambient 


Pa 


10~^gcm-3 


(*) sound speed of ambient 


Cs,a 


10~^ cms~^ 


(*) viscosity of ambient 




lO^^gcm"^ s~^ 
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COMPUTATIONAL DOMAIN 

- Cartesian coordinates (x,y,z) 

- Rest frame of mass center of droplet 




DROPLET 

- incompressible fluid 



AMBIENT REGION 

- Very low density 

- Compressible 



y, z 



RAREFIED GAS FLOW 

- Free molecular flow 

- Ram pressure on droplet surface 



Figure 1: Schematic picture of our numerical model and the coordinate system. 
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Figure 2: The staggered grid adopted in our model. The scalar variables are 
defined at the cell center and the velocity u = {u,v,w) is defined on the cell- 
edge. The external force like the surface tension and the ram pressure of the 
gas flow are defined at the cell-center. 
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^/'"o ^ 30 c ^/"^o — 30 cm 3 ' 

Figure 3: Time evolution of molten droplet exposed to the gas flow. The gas 
flow comes from the left side of panels. The initial droplet radius is rp — 500 /im 
(it corresponds to the Weber number We = 0.5). 
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(a) To = 100 Mm {W^ = 0.1) 



(d) To = 1000 /im (VK< 
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calculation 
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ave. = 0.966 - 



time [x 10' sec] 

(b) ro = 200 /xm (We = 0.2) 
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time !x 10'^ sec] 

(c) ro = 500 Atm {We = 0.5) 





5 10 15 20 25 30 35 40 45 
time txW' sec] 

(e) ro = 2000 /im (We = 2.0) 




10 20 30 40 50 60 70 
time fx 10'^ sec] 

(f) ro = 5000 Atm (We = 5.0) 




20 40 60 80 100 120 140 
time [x 10'^ seel 



Figure 4: Vibrational motions of molten droplets for various radii: the deforma- 
tion by the ram pressure and the recovery motion by the surface tension. The 
horizontal axis is the time since the ram pressure begins to affect the droplet and 
the vertical axis is the axial ratio of the droplet C/B (see text). The solid curves 
are the computational results and the dashed lines indicate the time-averaged 
values. 
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Sekiya et al. (2003) 

extrapolation 

ave. of calculation • 



0.8 
0.6 
0.4 
0.2 




~\ 1 — I — s~l — I I I 



Pfm = 4000 dyne cm' 
y = 400 dyne cni'^ 
=1.3 gem s 



J I I I I I I I 



_i I I I I I I I 



200 



2000 

droplet radius rg l\im] 



10000 



Figure 5: Time-averaged values of axial ratios C/B as a function of the initial 
droplet radius tq (filled circles). The linear solutions derived by Sekiya et al. 
(2003) are also displayed (solid curve). The dashed curve indicates a simple 
extrapolation of the linear solutions. 
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^/''o — 50 c ^/""o — 50 cm s ' 

Figure 6: Same as Fig. [3] except for ro — 5000 ^ni {We = 5). 
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(a) 0.068 sec 



(b) 0.143 sec 
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not fragment O 
cavitation • 
fragment ^ 




1 10 100 

coefficient of viscosity of droplet |Jd [g cm'^ s"^] 

Figure 8: Fates of the dust particles when they melt in the high- velocity gas 
flow. The horizontal axis is the coefficient of viscosity of droplet /id and the 
vertical axis is the droplet radius ro. We assume pfm = 4000 dyne cm^^, 7 = 
400 dyne cm^^, and pd — 3gcm~'^, respectively. Our simulation results are also 
shown (symbols). 
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U 




0.8 h 
0.75 
0.7 



Y = 400 dyne cm \ 



1 -1 



M-d = 1.3 gem s 



j_2 I I I I I I I \ I I I i. 



100 



1000 

droplet radius rg [\im] 



10000 



Figure 9: The axial ratios C/B of group- A chondrules are plotted as a function 
of the chondrule radius; filled circles are oblate (B/A > C/B) and open circles 
are others. The solid curves are the linear solutions for various ram pressures 
(Sckiya et al. 2003). The numbers in the panel indicate the values oipha in the 
unit of dynecm"^. The dashed curves are simple extrapolations of the linear 
solutions. 
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(a) low dust-to-gas mass ratio 



(b) high dust-to-gas mass ratio 




Figure 10: Dynamical/thermodynamical evolutions of the precursor dust par- 
ticles in the post-shock region are plotted as a function of the distance from 
the shock front. Top: the dust temperature. Td, and the radiation temperature 
Trad- Middle: velocities of dust particle, Vd, and gas Vg. Bottom: the ram 
pressure pfm- The time for dust particle after passage of the shock front is also 
displayed. 
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i-1 H ■^i+l 



Figure 11: Interpolate functions with various methods: CIP (sohd), Lax- 

Wendroff (dashed), and first-order upwind (dotted). The filled circles indicate 
the values of / defined on the digitized grid points Xi-\, Xi, and Xi+i before 
updated. 
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Figure 12: The schematic picture of the CIP-CSL2 scheme. 
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Figure 13: Numerical solutions of the one-dimensional advection or eonservative 
equation solved by various schemes: (a) first-order upwind, (b) Lax-Wendroff, 
(c) CIP, (d) R-CIP, (e) R-CIP-CSL2 without anti-diffusion, and (f ) R-CIP-CSL2 
with anti-diffusion. 
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